(*

**************************
* WRITTEN BY SIJA - 2024 *
**************************

********************************************************************************
* SOME EXAMPLES OF PROCEDURES DEFINITIONS FROM CLASSICAL DIFFERENTIAL GEOMETRY *
********************************************************************************

**********
* Curves *
**********


01. Tan_vec3D := proc(X, unon)
02. Nor_vec3D := proc(X, unon)
03. Bin_vec3D := proc(X, unon)

04. curvature_3D := proc(X)
05. curvatureR_3D := proc(X)

06. torshion_3D := proc(X)
07. torshionR_3D := proc(X)

08. curveTanPln_3D := proc(X, u, w)
09. curveTanPlnVec_3D := proc(X, tanVec, norVec, u, w)

10. curveNorPln_3D := proc(X, u, w)
11. curveNorPlnVec_3D := proc(X, binVec, norVec, u, w)

12. curveRecPln_3D := proc(X, u, w)
13. curveRecPlnVec_3D := proc(X, binVec, tanVec, u, w)

14. sphereRadius := proc(X)
15. sphereCenter := proc(X)
16. osculSphere := proc(X)

17. osculCircleRadius := proc(X)
18. osculCircleCenter := proc(X)

19. Circle_3D := proc(circleCenter, circleRadius, v1, v2)

20. rollCircCentInTanPln_3D := proc(X, r, inORout)
21. rollPntOnCircInTanPln_3D := proc(X, r, inORout)

18. rollCircCentInRecPln_3D := proc(X, r, inORout)
19. rollPntOnCircInRecPln_3D := proc(X, r, inORout)



************
* Surfaces *
************



*)


#**********************
# Curves - Procedures *
#**********************


Tan_vec3D := proc(X, unon) # X <- Column Vector
  local XL, var, Xd, XLx_d, XLy_d, XLz_d; 
      
  XL := convert(X, list);
  var := op(1, indets(X,name));
    
  XLx_d := diff(XL[1], var);
  XLy_d := diff(XL[2], var);
  XLz_d := diff(XL[3], var);
  
  Xd := < XLx_d, XLy_d, XLz_d >;
  
  # Column Vector of Tan_vec3D
  if unon
  then 
    simplify(eval(Xd/norm(Xd, 2))) assuming var::real;
  else   
    simplify(eval(Xd)) assuming var::real;
  end if;       
end proc:


Nor_vec3D := proc(X, unon) # X <- Column Vector
  local XTV, XL, var, Xdd, XLx_dd, XLy_dd, XLz_dd, vcp, vcpf; 
  
  XTV := Tan_vec3D(X, unon);
      
  XL := convert(X, list);
  var := op(1, indets(X,name));
      
  XLx_dd := diff(XL[1], var, var);
  XLy_dd := diff(XL[2], var, var);
  XLz_dd := diff(XL[3], var, var);
  
  Xdd := < XLx_dd, XLy_dd, XLz_dd >;
  
  vcp := VectorCalculus:-CrossProduct(XTV, Xdd);
  vcpf := VectorCalculus:-CrossProduct(vcp, XTV);
  
  # Column Vector of Nor_vec3D
  if unon
  then
    simplify(eval(vcpf/norm(vcpf, 2))) assuming var::real;
  else   
    convert(simplify(eval(vcpf)), Vector) assuming var::real;
  end if;       
end proc:


Bin_vec3D := proc(X, unon) # X <- Column Vector
  local XTV, XNV, XL, var, Xdd, XLx_dd, XLy_dd, XLz_dd, vcp, vcpf; 
  
  var := op(1, indets(X,name));
  
  XTV := Tan_vec3D(X, unon);
  XNV := Nor_vec3D(X, unon);
   
  vcp := VectorCalculus:-CrossProduct(XTV, XNV);
    
  # Column Vector of Bin_vec3D
  if unon
  then   
    simplify(eval(vcp/norm(vcp, 2))) assuming var::real;
  else
    convert(simplify(eval(vcp)), Vector) assuming var::real;
  end if;       
end proc:


curvature_3D := proc(X) # X <- Column Vector
  local XTV, XL, var, XLx_dd, XLy_dd, XLz_dd, Xdd, vcp; 
       
  XTV := Tan_vec3D(X, false);   
  
  XL := convert(X, list);
  var := op(1, indets(X,name));
      
  XLx_dd := diff(XL[1], var, var);
  XLy_dd := diff(XL[2], var, var);
  XLz_dd := diff(XL[3], var, var);
  
  Xdd := < XLx_dd, XLy_dd, XLz_dd >;
  
  vcp := VectorCalculus:-CrossProduct(XTV, Xdd);
  
  # Curvature in 3D of vector function X
  simplify(eval(norm(vcp, 2)/norm(XTV, 2)^3)) assuming var::real;
     
end proc:


curvatureR_3D := proc(X) # X <- Column Vector
  local var, cur3d; 
    
  var := op(1, indets(X,name));
  
  cur3d := curvature_3D(X);
    
  # Curvature Radius in 3D of vector function X
  if cur3d <>0 
  then 
    abs(simplify(eval(1/cur3d))) assuming var::real; 
  else 
    print("The radius of curvature is infinitely large");
    infinity; 
  end if;
     
end proc:


torshion_3D := proc(X) # X <- Column Vector
  local XL, var, Xd, XLx_d, XLy_d, XLz_d,
                 Xdd, XLx_dd, XLy_dd, XLz_dd,
                 Xddd, XLx_ddd, XLy_ddd, XLz_ddd,
                 vcp, dop, vcpd; 
      
  XL := convert(X, list);
  var := op(1, indets(X,name));
    
  XLx_d := diff(XL[1], var);
  XLy_d := diff(XL[2], var);
  XLz_d := diff(XL[3], var);
  
  Xd := < XLx_d, XLy_d, XLz_d >;
  
  XLx_dd := diff(XLx_d, var);
  XLy_dd := diff(XLy_d, var);
  XLz_dd := diff(XLz_d, var);
  
  Xdd := < XLx_dd, XLy_dd, XLz_dd >;
  
  XLx_ddd := diff(XLx_dd, var);
  XLy_ddd := diff(XLy_dd, var);
  XLz_ddd := diff(XLz_dd, var);
  
  Xddd := < XLx_ddd, XLy_ddd, XLz_ddd >;
  
  vcp := VectorCalculus:-CrossProduct(Xdd, Xddd);
  dop := VectorCalculus:-DotProduct(Xd, vcp);
  vcpd := VectorCalculus:-CrossProduct(Xd, Xdd);
  
  # Torshion in 3D of vector function X
  simplify(eval(dop/norm(vcpd, 2)^2)) assuming var::real;
  
end proc:
 

torshionR_3D := proc(X) # X <- Column Vector
  local var, tor3d; 
    
  var := op(1, indets(X,name));
  
  tor3d := torshion_3D(X); 
    
  # Torshion Radius in 3D of vector function X
  if tor3d <>0 
  then 
    abs(simplify(eval(1/tor3d))) assuming var::real; 
  else 
    print("The radius of torsion is infinitely large");
    infinity; 
  end if;
  
end proc:


curveTanPln_3D := proc(X, u, w) # X <- Column Vector
  local var, tanvec, norvec;
  
  var := op(1, indets(X,name));
      
  tanvec:= Tan_vec3D(X, true);
  norvec:= Nor_vec3D(X, true);
      
  # Parametric tangent plane to X 
  simplify(eval([X[1] + u*tanvec[1] + w*norvec[1], 
                 X[2] + u*tanvec[2] + w*norvec[2], 
                 X[3] + u*tanvec[3] + w*norvec[3]])) assuming var::real;
                
end proc:


curveTanPlnVec_3D := proc(X, tanVec, norVec, u, w) # X, tanVec, norVec <- Column Vectors functions of the same one variable
  local var;
  
  var := op(1, indets(X,name));
     
  # Parametric tangent plane to X 
  simplify(eval([X[1] + u*tanVec[1] + w*norVec[1], 
                 X[2] + u*tanVec[2] + w*norVec[2], 
                 X[3] + u*tanVec[3] + w*norVec[3]])) assuming var::real;
                
end proc:


curveNorPln_3D := proc(X, u, w) # X <- Column Vector
  local var, binvec, norvec;
  
  var := op(1, indets(X,name));
      
  binvec:= Bin_vec3D(X, true);
  norvec:= Nor_vec3D(X, true);
      
  # Parametric normal plane to X 
  simplify(eval([X[1] + u*binvec[1] + w*norvec[1], 
                 X[2] + u*binvec[2] + w*norvec[2], 
                 X[3] + u*binvec[3] + w*norvec[3]])) assuming var::real;
                
end proc:


curveNorPlnVec_3D := proc(X, binVec, norVec, u, w) # X, tanVec, norVec <- Column Vectors functions of the same one variable
  local var;
  
  var := op(1, indets(X,name));
     
  # Parametric normal plane to X 
  simplify(eval([X[1] + u*binVec[1] + w*norVec[1], 
                 X[2] + u*binVec[2] + w*norVec[2], 
                 X[3] + u*binVec[3] + w*norVec[3]])) assuming var::real;
                
end proc:


curveRecPln_3D := proc(X, u, w)  # X <- Column Vector
  local var, binvec, tanvec;
  
  var := op(1, indets(X,name));
      
  binvec:= Bin_vec3D(X, true);
  tanvec:= Tan_vec3D(X, true);
      
  # Parametric rectifying plane to X 
  simplify(eval([X[1] + u*binvec[1] + w*tanvec[1], 
                 X[2] + u*binvec[2] + w*tanvec[2], 
                 X[3] + u*binvec[3] + w*tanvec[3]])) assuming var::real;
                
end proc:


curveRecPlnVec_3D := proc(X, binVec, tanVec, u, w) # X, tanVec, norVec <- Column Vectors functions of the same one variable
  local var, binvec, norvec;
  
  var := op(1, indets(X,name));
 
  # Parametric rectifying plane to X 
  simplify(eval([X[1] + u*binVec[1] + w*tanVec[1], 
                 X[2] + u*binVec[2] + w*tanVec[2], 
                 X[3] + u*binVec[3] + w*tanVec[3]])) assuming var::real;
                
end proc:


sphereRadius := proc(X) # X <- Column Vector
  local XL, var, XLx_d, XLy_d, XLz_d, Xd, cur3D, cur3D_d, 
        tor3d, c_b; 
    
  XL := convert(X, list);
  var := op(1, indets(X,name));
    
  XLx_d := diff(XL[1], var);
  XLy_d := diff(XL[2], var);
  XLz_d := diff(XL[3], var);
  
  Xd := < XLx_d, XLy_d, XLz_d >;
           
  cur3D := curvature_3D(X);
  cur3D_d := diff(cur3D, var);
  tor3d := torshion_3D(X);
  
  simplify(eval(sqrt(1/cur3D^2 + cur3D_d^2/(norm(Xd, 2)^2*cur3D^4*tor3d^2)))) assuming var::real;
 
end proc:


sphereCenter := proc(X) # X <- Column Vector
  local XL, var, XLx_d, XLy_d, XLz_d, Xd, 
        cur3D, cur3D_d, tor3d, XBV, XNV, c_b;
        
  XL := convert(X, list);
  var := op(1, indets(X,name));
    
  XLx_d := diff(XL[1], var);
  XLy_d := diff(XL[2], var);
  XLz_d := diff(XL[3], var);
  
  Xd := < XLx_d, XLy_d, XLz_d >;
  
  cur3D := curvature_3D(X);
  cur3D_d := diff(cur3D, var);
  
  tor3d := torshion_3D(X);
  
  XBV := Bin_vec3D(X, true);
  XNV := Nor_vec3D(X, true);
  
  # Osculation sphere center
  simplify(eval(X + (1/cur3D)*XNV - cur3D_d/(norm(Xd, 2)*cur3D^2*tor3d)*XBV)) assuming var::real;
 
end proc:


osculSphere := proc(X) # X <- Column Vector
  local var, sC, sR, sp;
  global theta, phi;

  var := op(1, indets(X,name));  

  sC := sphereCenter(X); 
  sR := sphereRadius(X);
  # Attention: theta and phi are globa variables
  sp := convert([cos(theta)*sin(phi)*sR, sin(theta)*sin(phi)*sR, cos(phi)*sR], Vector)   assuming var::real;
    
  # Osculation sphere           
  simplify(eval(sC + sp));
  
end proc:


osculCircleRadius := proc(X) # X <- Column Vector
  local var, cur3d; 
    
  var := op(1, indets(X, name));
  
  cur3d := curvature_3D(X);
    
  # Osculation Circle Radius in 3D of vector function X
  if cur3d <>0 
  then 
    abs(simplify(eval(1/cur3d))) assuming var::real; 
  else 
    print("The radius of the osculation circle is infinitely large");
    infinity; 
  end if;

end proc:


osculCircleCenter := proc(X) # X <- Column Vector
  local var, cR, unvec;
  
  var := op(1, indets(X,name));
  
  cR:= osculCircleRadius(X);
  unvec:= Nor_vec3D(X, true); 
 
  # Osculation circle center          
  simplify(eval(X + cR*unvec)) assuming var::real;
     
end proc:


Circle_3D := proc(circleCenter, circleRadius, v1, v2)
  local var, xc, yc, zc, center, radius;
  
  var := op(1, indets(circleCenter,name));
  
  center := circleCenter; 
  radius := circleRadius;
    
  # Parametrize circle
  # s := 0..2*Pi is a circle parametr:
  xc := center[1] + radius*cos(s)*v1[1] + radius*sin(s)*v2[1];
  yc := center[2] + radius*cos(s)*v1[2] + radius*sin(s)*v2[2];
  zc := center[3] + radius*cos(s)*v1[3] + radius*sin(s)*v2[3];
  
  # Cirrcle on a plane spread over vectors v1 and v2 in 3D
  simplify(eval([xc, yc, zc])) assuming var::real;
  
end proc:


rollCircCentInTanPln_3D := proc(X, r, inORout) # X <- Column Vector
  local var, unvec;
  
  var := op(1, indets(X, name));
  
  unvec:= Nor_vec3D(X, true);
   
  if inORout 
  then
    # A circle with this center rolls along the "concavity" - in of the curve
    simplify(eval(X + r*unvec)) assuming var::real;
  else
    # A circle with this center rolls on the "convexity" - out of the curve   
    simplify(eval(X - r*unvec)) assuming var::real;
  end if;
  
end proc:


rollPntOnCircInTanPln_3D := proc(X, r, inORout) # X <- Column Vector
  local XL, XLx_d, XLy_d, XLz_d, Xd, var, tanvec, norvec, rollCC, 
        arcelem, rarc, rotang;
  
  XL := convert(X, list);
  var := op(1, indets(X,name));
    
  XLx_d := diff(XL[1], var);
  XLy_d := diff(XL[2], var);
  XLz_d := diff(XL[3], var);
  
  Xd := < XLx_d, XLy_d, XLz_d >;
  
  arcelem := sqrt(VectorCalculus:-DotProduct(Xd, Xd));
  arcelem := subs(var = u, arcelem);
  
  rotang:= evalf(Int(arcelem, u = 0..var)/r) assuming var::real;
      
  tanvec := Tan_vec3D(X, true);
  norvec := Nor_vec3D(X, true);
  rollCC := rollCircCentInTanPln_3D(X, r, inORout);
              
  # Rolling Point On Circle In Tangent Plain in 3d
  if inORout
  then   
    rollCC - r*tanvec*cos(rotang) + r*norvec*sin(rotang);
  else
    rollCC + r*tanvec*sin(rotang) - r*norvec*cos(rotang);
  end if;

end proc:


rollCircCentInRecPln_3D := proc(X, r, inORout) # X <- Column Vector
  local var, ubvec;
  
  var := op(1, indets(X, name));
  
  ubvec:= Bin_vec3D(X, true);;
  
  if inORout 
  then
    # A circle with this center rolls along the "concavity" - in of the curve
    simplify(eval(X + r*ubvec)) assuming var::real;
  else
    # A circle with this center rolls on the "convexity" - out of the curve   
    simplify(eval(X - r*ubvec)) assuming var::real;
  end if;
  
end proc:


rollPntOnCircInRecPln_3D := proc(X, r, inORout) # X <- Column Vector
  local XL, XLx_d, XLy_d, XLz_d, Xd, var, tanvec, binvec, rollCC, 
        arcelem, rarc, rotang;
  
  XL := convert(X, list);
  var := op(1, indets(X,name));
    
  XLx_d := diff(XL[1], var);
  XLy_d := diff(XL[2], var);
  XLz_d := diff(XL[3], var);
  
  Xd := < XLx_d, XLy_d, XLz_d >;
  
  arcelem := sqrt(VectorCalculus:-DotProduct(Xd, Xd));
  arcelem := subs(var = u, arcelem);
  
  rotang:= evalf(Int(arcelem, u = 0..var)/r) assuming var::real;
    
  tanvec := Tan_vec3D(X, true);
  binvec := Bin_vec3D(X, true);
  rollCC := rollCircCentInRecPln_3D(X, r, inORout);
            
  # Rolling Point On Circle In Tangent Plain in 3d
  if inORout
  then   
    rollCC - r*tanvec*cos(rotang) + r*binvec*sin(rotang);
  else
    rollCC + r*tanvec*sin(rotang) - r*binvec*cos(rotang);
  end if;
  
end proc:
